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Abstract 

The aim of this paper is to compare resuhs from lattice-Boltzmann and Brownian dynamics 
simulations of linear chain molecules. We have systematically varied the parameters that may 
affect the accuracy of the lattice-Boltzmann simulations, including grid resolution, temperature, 
polymer mass, and fluid viscosity. The effects of the periodic boundary conditions are minimized by 
an analytic correction for the different long-range interactions in periodic and unbounded systems. 
Lattice-Boltzmann results for the diffusion coefficient and Rouse mode relaxation times were found 
to be insensitive to temperature, which suggests that effects of hydrodynamic retardation are 
small. By increasing the resolution of the lattice-Boltzmann grid with respect to the polymer size, 
convergent results for the diffusion coefficient and relaxation times were obtained; these results 
agree with Brownian dynamics to within 1-2%. 
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I. INTRODUCTION 



The fluctuating lattice-Boltzmann (FLB) equation has been proposed as a basis for 



numerical simulations of polymer solutions 



2|, using a frictional coupling between the poly- 



mer and the surrounding fluid. Single molecule dynamics obtained with this algorithm 
compared favorably with coarse-grained molecular dynamics simulations using explicit sol- 
vent [3]. More recently, the same idea has been investigated for confined polymers |[^], as 
a simpler and possibly more efficient alternative to Brownian dynamics (BD) Both 
methods have been applied to problems of polymer migration in shear and pressure-driven 
flows and showed similar trends with increasing shear rate, but different choices of 

polymer model precluded a quantitative comparison. Here we initiate a systematic compar- 
ison of FLB and BD methods, beginning with the properties of an isolated chain; subse- 
quently we will extend the investigation to confined polymers in shear and pressure-driven 
fiows. This work complements a recent study by Pham et al. [sj]. 

Inertial effects are neglected in Brownian dynamics, resulting in instantaneous propa- 
gation of momentum. Although the solvent degrees of freedom are thereby eliminated, 
the interactions between the beads are long-range, which leads to an 0{N^) scaling of the 
computational cost for a polymer consisting of segments. In the FLB method, all in- 
teractions are local and the computational cost scales linearly with the volume. However, 
the lattice-Boltzmann method introduces an extra, inertial time scale, during which the 
hydrodynamic interactions propagate throughout the fluid by viscous momentum diffusion. 
Surprisingly this has little effect on the time step; FLB simulations use comparable time 
steps to BD, as will be seen in Sec. IIII A[ Hydrodynamic retardation, sometimes thought 
to be a potential source of error in lattice-Boltzmann simulations of Stokes flow 0, ^ , is in 
fact easily managed. Nevertheless, the additional degrees of freedom of an explicit solvent 
model adds considerably to the computational cost. Generally, dilute systems in unconflned 
geometries favor BD, while more concentrated solutions in conflned geometries favor the 
FLB method 

Quantitative comparisons require an identical micro-mechanical model of the polymer. 
However, the BD simulations are for an isolated chain, while the FLB simulations use 
periodic boundary conditions. It is therefore necessary to correct for the differences in 
the long-range flow flelds in periodic and unbounded systems; in particular, the diffusion 
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coefficient in a periodic system has a correction proportional to L ^, where L is the length 



of the periodic unit cell. However, prior research [4, [lOl] has shown that these corrections can 
be calculated quantitatively, based on the hydrodynamic theory for a periodic unit cell 
Corrections to the configurational properties and Rouse relaxation times are much smaller, 
0{L^^)\ when the box length is 5 — 10 times the polymer size, deviations from the infinite 
system are negligible. The diffusion coefficients and Rouse relaxation times were found to 
depend on the degree of discretization of the lattice-Boltzmann fluid. Our results show 
numerical convergence with increasing grid resolution, and the converged results agree with 
Brownian dynamics within 1 — 2%. 



II. POLYMER MODEL AND SIMULATION METHODS 

The polymer model consists of + 1 beads connected by FENE springs between neigh- 
boring beads, 

A 1 / r2\ 

$5 = X.'^sdrm -^il), (t)s{r) = --nrl\Yi[l- — ], (1) 

where k is the spring constant, tq is the maximum extension of the spring, and fi is the 
position vector of the i*^ bead. The simulations use a value of tq = 5.486, where h = a/T/k 
and T is the thermal energy corresponding to an absolute temperature T/kB- The root- 
mean-square bond length in an ideal chain, (r^)^^^ = 1.606, follows from the potential in 

Eq. O; 



In addition to the FENE potential, there is an excluded volume interaction between the 
beads, 

^EV = ^(pEvi\ri-rj\), 0ijy(r) = eexp(-/?r^), (3) 
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with p = 1.506"^ and e = 2. TIT. The potential ener gy $ = + ^ev approximates a 
DNA molecule with ~ 10 Kuhn segments per spring [sl-ll^. The identical polymer model, 
described by Eqs. ([T]) and ([3]), was used for both BD and FLB simulations. The beads are 
coupled to the fluid with a Stokes friction coefficient ^ = Qnrja, where rj is the fluid viscosity 
and the hydrodynamic radius of the beads a = 0.3626. 
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A. Brownian dynamics 



Brownian dynamics neglects inertia, and the state of the polymer is therefore completely 
specified by the positions of the + 1 beads, rj. Hydrodynamic interactions (HI) between 
the beads are introduced in a pairwise-additive approximation through the mobility matrix 
/Xjj, which connects the mean (or drift) velocity of bead i to the force on bead j, 



N 



Vi 



(4) 



The conservative force, Fa 



j=0 

Vt- $, is derived from the potential energy of the polymer. 



Eqs. (PP) and (El). We use the Rotne-Prager regularization of the mobility matrix of point 



particles 13|, [ij], which approximates the Stokes-flow result in the limit Vij ^ a, but ensures 



that the mobility matrix remains positive definite for all rij: 



C[I + 



> 2a 



— , < 2a 



(5) 



where 





3 a la^ 


C2 




4 rij 2 rf- ' 




C[ = 


q J,.. 
1-— ^, 




32 a 





3 a 
32 a ' 



3 a^ 
2^' 



(6) 



When i = j, the self mobility, = ^ ^I, is the mobility of an isolated sphere. 

The Rotne-Prager mobility is divergence free, ' ^^ij ~ ^^^^ order 

[ ^ 

Ermak and McCammon algorithm \15^ for stochastic integration reduces to an explicit Euler 
integration scheme, 

ri{t + At) = Viit) + ViAt + Awi, (7) 
where At is the time step and Awi is a random displacement with zero mean and covariance 

< Awi (g) Awj >= 2T/x.^.At. (8) 
We used a Cholesky decomposition of the grand mobility matrix. 



t^Nfi Mat,! 



l-''N,N 



(9) 



4 



to calculate the random displacements, which is an 0{N^) computation. However, the 
Cholesky decomposition is a small fraction of the total computational cost for the short 
chains {N = 10, 20) used in this work, and eliminates any possible errors associated with 
the Chebyshev polynomial approximation [16], which scales more favorably as 0(A^^'^^). 



B. Lattice Boltzmann 



The fluctuating lattice-Boltzmann model 



of dilute polymer solutions in periodic 
original formulation of the FLB model 




has been used to simulate the dynamics 
19| and confined geometries ^, Q, ■ In the 
the viscous stress tensor was assumed to 



fluctuate around the local Navier-Stokes stress, but this model fails to satisfy the fluctuation- 
dissipation theorem at small scales |2l| unless thermal fluctuations in the non-hydrodynamic 
modes are included as well. Here we summarize the improved method, following the recent 
statistical-mechanical formulation of the FLB equation 22]. Including thermal fluctuations 
in the non-hydrodynamic modes leads to small, but noticeable, improvements in the equipar- 
tition of energy between the fluid and polymer degrees of freedom (Sec. IIII A|) . 

In the lattice-Boltzmann (LB) model, the fluid degrees of freedom are represented by 
a discretized one-particle velocity distribution function ni{r,t), which describes the mass 
density of particles with velocity Cj at the position r and time t. The hydrodynamic fields, 
mass density p, and momentum density j = pu, are moments of this velocity distribution. 

The time evolution of ni{r,t) is described by a discrete analogue of the Boltzmann equa- 
tion [23], 

riiir + CiAt, t + At)= ni{r, t) + [n(r, t)] ; (11) 

here Aj is the change in due to instantaneous collisions at the lattice nodes and At is the 
time step. 

The D3Q19 model 2J] was used, which includes rest particles and 18 velocities corre- 
sponding to the [100] and [110] directions of a simple cubic lattice. The population density 
associated with each velocity has a weight a"^' that describes the fraction of particles with 
velocity Cj in a fluid at rest: 



1 

3' 



1 



1 

36' 



(12) 



The deterministic collision operator is typically linearized about the low-velocity equilibrium 
distribution, n^*^, 



(p, u) = a^^p 1 + 



+ 

where the speed of sound = 3~^/^ Ax/At and Ax is the lattice spacing. T 



(13) 



ne non- 



equilibrium distribution, n^^'' = nt — n^'^, can then be expanded in moments |25l. l26l| . 



TTlk 



(14) 



using tensorial polynomials of the lattice vectors, basis. We use a different basis 



from Refs. 



25 



26| . such that the back transformation includes the weights a^^ 27l |. 

'^w^^Ckimk, (15) 



where = 'Yli^'^^^'ki normalizing factor for e^,. During the collision process, the 

moments relax towards equilibrium (zero). 



(16) 



where the relaxation parameter is bounded by 1 7^1 < 1. In these simulations we used a 
two-parameter collision operator, with different eigenvalues for the modes with odd (70) and 
even (7e) powers of Ci [28|. The shear viscosity is related to 7e, 

pc^h / 1 + 7e 



and 



7o 



l-7e 
77e + l 



(17) 



(18) 



271, 



28|, 



7e + 7' 

this relation makes the location of a planar solid boundary independent of viscosity 
although that property is not essential in the present context. 

The key difference between the FLB and LB models is in the collision operator. The 
FLB collision operator contains random excitations of the non-conserved moments, c.f. 
Eq. ([T6D y, 

/ rvm 11)1 (A — o/?^ 

(19) 



pmpWk{l-ll) 
k - Ikrrik + \l -Vk^ 



Ax3 



where (pk is a random variable with zero mean and unit variance. It is important to use a 



bounded distribution of random numbers 



29|, or large changes in nik will occasionally occur, 



leading to negative values of n^. The amplitude of the random forcing is determined from 



18| and is controlled by the mass of an LB "particle", 



the fluctuation-dissipation relation 
rup [22] . Thermodynamic consistency requires that rrip is related to the effective temperature 



of the fluctuating fluid 



22 



rripcl 



T. 



(20) 



In this work the random forcing is applied ^all the non-conserved modes 2l|, not just the 
stress jl]. It has been shown theoretically [22] that the excitation of the non-hydrodynamic 
modes is essential to satisfy the fluctuation-dissipation relation at all scales, although at 
long wavelengths only the excitations in stress are important. 



C. Polymer-fluid coupling 

The polymer is coupled to the LB fluid by frictional forces between the beads and the 
fluid. The equations of motion for the i*^ monomer can be written in inertial form as, 

^ = Vi, = Fi-^o [Vi{t) - uin, t)] + R,{t). (21) 

at at 

where the hydrodynamic force includes a frictional drag, based on the difference in velocity 
between the bead and the surrounding fluid, and a random force, Ri, to balance the addi- 
tional dissipation [2]. Since the fluid satisfles its own fluctuation-dissipation relation, Ri has 
a local covariance matrix 

{R,{t)R,{t')) = 2T^oS{t - t')6,,I. (22) 

Hydrodynamic interactions between the beads are transmitted through the fluid via corre- 
lated fluctuations in the velocity fleld, which develop over the inertial time scale, pr'^/f], where 
r is the separation between beads. The large time-scale separation between the dynamics 
of the polymer and the individual monomers allows time for the hydrodynamic interactions 
to reach a quasi-steady state, without imposing this condition at each and every time step. 
We will show numerically that both inertial (FLB) and diffusive (BD) simulations can use 
similar time steps, of the order of the monomer diffusion time (see also Ref. 

Since the monomers move continuously over the grid, the fluid velocities, it„, are inter- 
polated from neighboring grid points to the bead location r^, 

tt(ri,t) = ^ A(ri - r„)it„. (23) 

n 
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The interpolating function A{rx,ry,rz) is taken as a product of one-dimensional func- 



tions 



30| 



A{x,y,z) = (f) 



(24) 



(25) 



where Ax is the grid spacing in the lattice-Boltzmann simulations. Typically the weights 
0(m) are determined by linear (two-point) interpolation, 

1 — \u\ \u\ < 1, 
|m| > 1, 

but a more precise interpolation is possible using three or four points in each coordinate 
direction. Numerical tests of the different interpolations can be found in Ref. 3|. We will 
present some results with three-point interpolation, 



I (1 + VI - Su^) 
1(5- 3\u\ - v/-2 + 6|m| -3^2^ 







0<\u\< h 



< u < 



(26) 



but most of the results use 2-point (linear) interpolation. To conserve momentum, the 
accumulated force exerted by the bead on the fluid is distributed to the surrounding nodes 



with the same weight function 

The input friction ^0 = Svrr^ao is not the same as the effective friction C, = Qrcrja, as 

measured by the drag force on the bead or by its diffusion coefficient. This is because the 

force added to the fluid renormalizes the input friction, 

1 



where g is a. numerical factor 



1 _ 1 

? ^0 QnijAxg '' 



(27) 



29| that is independent of the fluid viscosity but depends on the 



interpolation function. From the diffusion of individual monomers we have determined values 
of g = 1.3 for linear (two-point) interpolation and g = 1.0 for the three-point interpolation. 
The FLB results in this paper are matched to BD simulations with the same effective radius, 
a. 

The coupled equations of motion for the particles and fluid are solved by operator split- 
ting [29:]; typically, the thermodynamic forces are integrated with a smaller time step than 
the hydrodynamic forces to maintain stability. The LB time interval At is decomposed into 
M steps of length h = At/M, where M is chosen to be sufficiently large that the conserva- 
tive forces are integrated accurately; typically M ~ 10 in our simulations. The algorithm 
used in this work is as follows: 
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1. At the beginning of the LB step, determine the fluid velocities at the grid points. 

2. Update the polymer positions and velocities over M sub-cycles. For each subcycle a 
modified Verlet algorithm is used to update the positions and velocities of the beads: 

(a) First stream the particle positions and velocities for half a time step, 

r-f ^ = rm + \<tl (28) 
= ^.W + ^i^f \ (29) 
where the conservative force F\ = — V (1)$ is evaluated from the coordinates at 

i 

the half time step, rf'\ 

(b) Use the updated positions, r^l\ to interpolate the fluid velocity to the bead 
locations, Eq. fl2^ . 

(c) Update the bead velocities for a full step h, using a midpoint approximation to 
the frictional drag force Eq. (1211) . 



,30, 

where (p^ is a vector of bounded random numbers with zero mean and unit vari- 
ance. 

(d) Redistribute the momentum transferred by particle-fluid coupling 
back to the fluid. 

(e) Stream the particle positions and velocities for the second half step, 

v,{t + h) = vf^ + ^Ff\ (32) 

r,{t + h)=rf^ + ^v,{t + h), (33) 

(34) 

The exact sequence of updates is important to preserve the second-order accuracy 
of the operator-splitting method. The algorithm reduces to the Verlet scheme when 
^0-0. 



3. Update LB populations to account for momentum transfer from particle- fluid coupling. 



4. Update FLB algorithm for one time step At. 

There are a number of nearly equivalent ways to break down the coupled dynamics of the 
particle-fluid system; the algorithm described above is the most accurate of the variations we 
have investigated, although the differences in long-time properties (conformation, diffusion. 
Rouse relaxation times) are generally small. The midpoint method is preferable to a first- 
order update, either explicit or implicit, since neither of these lead to exact thermalization 
of the kinetic energy of the particles. For force-free particles it is straightforward to show 
that [4 1 



with the plus sign following from the implicit update and the minus sign from the explicit 
update. By contrast, the midpoint method gives m (vf) = 3T exactly. There is a choice as 
to whether the momentum transferred to the fluid, Eq. (13T1) . affects the interpolated fluid 
velocity after each sub step (h) or only after each LB step (At). Numerical results show that 
the polymer temperature is closer to the fluid temperature (within 0.3%) if the interpolated 
velocity is updated every sub step (h). When the velocity is only updated at the LB steps 
(every At), the polymer temperature differs from the fluid temperature by about 3%. The 
results presented in Sec. [1111 have the interpolated fluid velocity updated every h. 

III. RESULTS 

The purpose of these simulations was to make precise comparisons of BD and FLB re- 
sults for an identical polymer model. We have compared static properties (radius of gyration 
and end-to-end distance) and dynamic properties (diffusion coefficient and Rouse relaxation 
times). The effects of time step have been investigated, and, in the case of the FLB simula- 
tions, the effects of grid resolution, temperature (fluctuation amplitude), fluid viscosity, and 
bead mass as well. To obtain statistically precise data, every FLB data point was calculated 
from an ensemble average over 160 different initial conditions. Each sample was equilibrated 
for a time of approximately 500to and data was collected for a further SOOOto, for a total of 
8 X lO^to- The time unit to = ^/k, and the Zimm time tz = Qn^Rg/T = 56.7to- The Brow- 
nian dynamics simulations of diffusion and Rouse relaxation times were run for 1.6 x lO^to 




(35) 
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FIG. 1: Conformational properties, Rg/b (closed symbols) and R^/b (open symbols), versus the 
dimensionless time step At/to- The FLB data used a grid resolution Ax = 1.296; conformational 
properties with other grid resolutions are statistically indistinguishable (see Table H]). 



and 8 x lO^to respectively. 



A. Static properties 

The mean square radius of gyration, R^, 

and end-to-end vector, R^, 

{Rl) = {ir^-r,r), (37) 

of the polymer are compared in Fig. [H The conformational properties from lattice- 
Boltzmann are indistinguishable from Brownian dynamics within the statistical errors 
(0.5%). Despite the extra inertial time scale, the FLB simulations use comparable timesteps 
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to Brownian dynamics; neither method shows statistically significant deviations in Rg and 
Re when the time step At is less than O.Olto- 

The accuracy of a BD simulation depends only on the time step, but results from FLB 
simulations may depend on a number of parameters: fiuctuation level or temperature (T), 
length of the periodic unit cell (L), grid resolution (Ax), fiuid viscosity (r^), and particle 
mass (m). Results for a range of values of these parameters are summarized in Tabled using 
a FENE chain (Sec. [Tll) of 10 segments. In FLB simulations, the timescale to = C,/n = C.b'^/T 
is controlled by the temperature, which sets the level of fiuctuations in the fluid, Eqs. (fT9|) - 
( |20|) . and particles, Eqs. (l2Ti)-(l22l). The temperature reflects the degree of coarse graining 



of the molecular fluid rather than the thermodynamic properties of the chain [22|, |29||; one 
dimensionless measure is the parameter a = (uf) /cj. (see Table H]), which relates the fluid 
velocity fluctuations in a single LB cell to the sound speed. For the LB model to adequately 
represent an incompressible fluid {\u\/cs < 0.1), a should be less than 0.01. At the highest 
temperature shown in Table [H a = 0.024, the polymer is indeed slightly swollen. 

The fluctuation-dissipation theorem (FDT) should ensure that the polymer thermalizes 
to the same temperature as the fluid; in other words m (vf) = M (ul), where M = pAx^ is 
the mass of fluid in a single LB grid cell. The column AT/T in Table [T] measures the relative 
deviations in the polymer temperature from thermal equilibrium; these are usually small, of 
the order of 0.2%. Larger deviations occur when the temperature is too high {a = 0.024), or 
when the LB model is not exactly thermalized (footnotes 5 and 6). If the kinetic (or ghost) 
modes are not subject to random forcing, the fluctuation-dissipation relation is broken at 
short length scales 
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22| . This causes a small error in the size of the polymer, 1 



when compared to the properly thermalized simulations, with similar deviations in the 
diffusion coefficient and Rouse relaxation times (footnote 5). Larger errors in both static and 
dynamic properties occur if the fluid dynamics is purely dissipative (footnote 6), because the 
fluctuation-dissipation relation is then broken at all length scales; results from simulations 
without fluid fluctuations, for example Ref. 9|, are invalid. The establishment of good 
thermal equilibrium between the polymer and fluid requires exact thermalization of the LB 
fluid and the coupling algorithm described in Sec. Ill C[ 

The ratios k/T and e/T must be kept constant if the polymer conformations are to be 
independent of the degree of coarse graining of the fluid degrees of freedom. The dimension- 
less timestep At /to of the FLB simulation in Fig. [T]is then controlled by the temperature 
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44.3 ± 0.2 7.51 ± 0.01 


0.208 


oo 


16.5 



TABLE I: Static and dynamic properties of a polymer chain. Fluctuating LB simulations for a 10- 
segment chain are compared with Brownian dynamics. The resolution of the LB grid, Ax, can be 

1 /9 

compared with the RMS distance between the beads ^r^^ = 1.606. The parameter a = {uf^ /c^ 
is a measure of the temperature of the fluid, T = M (uf ), where M = pAx^ is the mass of fluid in 
a single grid cell. The dimensionless time step in the FLB simulations is related to the temperature 
through the scaling with to = ^/ k = £J}^ /T. AT is the difference in temperatures of the particles 
and fluid. The mass of a bead, kinematic viscosity of the fluid and the length of the periodic unit 
cell are m = O.IM, v = O.lAx^/At and L = 9.2Rg, unless otherwise indicated. The statistical 
errors in diffusivity, D/Dq, and Rouse-mode relaxation time, ri/to, are less than 0.5%; Dq = T/^ 
is the monomer diffusivity, and Scq and Sc are the Schmidt numbers based on the monomer and 
polymer diffusivities. 

^ Length of periodic unit cell L = ISARg. 

2 The viscosity of the fluid is varied: (a) v = O.ObAx'^/At; (b) ly = 0.2Ax^/At. 
^ The mass of the bead is varied: (a) m = M; (b) m = lOM. 

^ Three-point interpolation. 

^ No excitation of the kinetic (ghost) modes. 

^ No excitation of the fluid modes. 
Results from BD simulations. 
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of the fluctuating fluid (T or a). Since the viscosity is independent of temperature in the 
FLB model, the Schmidt number Sc = rj/pD, varies inversely with T. The results in Ta- 
ble H] show that the static and dynamic properties are both insensitive to Schmidt number, 
but there are small deviations when Sc < 30, which is consistent with earlier findings {4]. 
Other authors have suggested that the Schmidt number based on the monomer diffusion, 
Scq = rj/pDo, should be in excess of 30 |2|], but our results suggest that this may be overly 
restrictive; we do not find systematic deviations in either static or dynamic properties until 
Scq < 10. Both static and dynamic properties are statistically independent of fiuid viscosity 
(footnote 2) and particle mass (footnote 3) over the ranges studied. Finally, we note that 
the conformational properties and Rouse-mode relaxation times are independent of the size 

n 

of the unit cell (footnote 1). These results are consistent with recent FLB simulations |8|, 
which found a weak system size dependence when L < 5Rg; in our simulations L ~ lORg. 
The systematic dependence of the diffusion coefficient on L has been analytically corrected 
in Table [H as discussed in detail in Sec. IIII B[ 



B. Diffusion coefficient 



The diffusion coefficient of the polymer is determined from the time-dependent displace- 
ment, Arc(t) = rc(t) — Tc(0), of the center of mass vector, Vc = {N + l)"^X]ilo^«- 
calculate the diffusion coefficient from the time derivative, 

D{t) = ~{Ar^{t)-Ar^{t)), (38) 

rather than the slope, since the derivative asymptotes at much earlier times. The short-time 
diffusivity determined by Brownian dynamics is found from Eq. (!38|) in the limit t — ^ 0. It is 
equal to the Kirkwood diffusivity and only slightly different, by 1 — 2%, from the long-time 
diffusivity 3l|- The FLB simulations are inertial, and here limj_+o-D(t) = 0. Nevertheless, 
in both methods the diffusivity reaches its asymptotic value, D, over a time of the order of 
the Zimm time, tz = GtttjR'^/T. 

Our investigations show that the diffusion coefficient in an FLB simulation depends on 
only three parameters: the size of the periodic unit cell, L, the grid resolution, Ax/b, and 
the method of interpolation. Since the BD results are for an isolated polymer, the FLB 
data must be corrected for finite-size effects; here we use a well-established correction for 
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0.25 



FIG. 2: Diffusion coefficients from FLB simulations with different unit cell sizes, L. The coarsest 
grid resolution, Ax/6 = 2.58, was used for computational efficiency. 



the diffusion coefficient [10|, which effectively eliminates the dependence of the diffusion 
coefficient on the box size. The self- diffusion coefficient of a spherical particle of radius R 
in a periodic system with repeat length L can be written as lOj, 



(39) 



where is the diffusion constant in a system of length L and = T/6iTr]R is the 
diffusion coefficient of an isolated sphere. At large distances, the average flow fleld around 
the polymer is similar to the flow fleld around a spherical particle; we therefore expect the 
same relation for the center-of-mass diffusion coefficient of the polymer. Although we do 
not know a priori what the effective radius of the polymer is, the leading order correction 
is independent of R, 



2.837T 
GirriL 



(40) 



The system-size dependence has been investigated using the coarsest resolution of the 
LB grid, Ax/b = 2.58, to maximize computational efficiency. The results in Fig. [2] show 
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L/Rg 




a 








2.8 





0015 


0.083 





216 


0.205 


4.7 





003 


0.114 





193 


0.191 


4.7 





0015 


0.115 





194 


0.192 


9.4 





003 


0.149 





189 


0.188 


9.4 





0015 


0.148 





188 


0.187 


18.9 





003 


0.168 





188 


0.188 


18.9 





0015 


0.168 





188 


0.187 


28.3 





003 


0.173 





186 


0.186 


28.3 





0015 


0.170 





184 


0.184 


28.3 





00075 


0.174 





188 


0.188 


28.3 





000375 


0.174 





188 


0.188 


37.7 





003 


0.173 





183 


0.183 


37.7 





0015 


0.175 





185 


0.185 


37.7 





00075 


0.176 





186 


0.186 


37.7 





000375 


0.177 





187 


0.187 



TABLE II: Effect of system size on the diffusion coefficient of a 10-segment chain. The resolution 
of the LB grid, l^x/h = 2.58. The parameter a = (uf) /c^ is a measure of the temperature of the 
fluid, T = M (uf ), where M = pAx^ is the mass of fluid in a single grid cell. Dl is the diffusion 
coefficient from FLB simulations and Doo is the corrected value from Eq (jlOp ; results including the 
(R/L)^ correction are indicated by 

the expected linear dependence on L~^, with the same asymptotic value of the polymer 
diffusivity {D/Dq = 0.1875) from either extrapolation, fitting to 4 different system sizes 
(9.5 < L/Rg < 38), or from a single simulation with L ~ lO-R^. Since the computational 
cost scales as L^, a single simulation takes 1/36 the time of a sequence of three simulations 
with box lengths in the ratio 1:2:3. Moreover, larger systems require a lower temperature 
to obtain the asymptotic (with T) diffusion coefficient, as shown in Table [ITl Data for large 
systems (L > 20Rg) shows that the limiting, low— T diffusion coefficient, as plotted in Fig. [2l 
requires a temperature 4 — 8 times smaller, which translates to 4 — 8 times more processing 
for the same statistics. 

It is possible to further refine the correction for finite-size effects by including the next 
term in Eq. fl39p . but this requires the polymer size. Defining the diffusivity of the isolated 
chain, Dqo = T/GnrjRao, in terms of the effective radius R^o, and rearranging Eq. fl5^ results 
in a cubic equation for x = Roo/L, 

4 

-7ra;3 - (2.837 + xl)x + 1 = 0, (41) 
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1 1.5 
Ax/b 



2.5 



FIG. 3: Diffusion coefficients from FLB simulations with different grid resolutions, Ax/b, are com- 
pared with Brownian dynamics (dashed line); the vertical bar indicates the statistical uncertainty 
in the BD data. 

where xl = Rl/L, with Dl = T/GnrjRL. Since the additional correction is small, Eq. (HTl) 
can be solved for a; in a few iterations. This leads to slightly more consistent diffusion 
coefficients from the smaller box sizes, as shown in final column of Table HTl The diffusivities 
in Table [T] include the extra correction term. 

Polymer diffusion coefficients from FLB simulations depend on grid resolution, Aa;, and 
previous work 0, 4| suggests that a ratio of Ax/b ~ 1 — 2 is adequate for most purposes. 
Here the diffusion coefficient of 10 segment chains are shown in Fig. [3] for a range of different 
grid resolutions, 0.65 < Ax/b < 2.58. The BD result is shown as a dashed line for clarity 
in comparison. By decreasing the ratio Ax/b, the number of grid points over which the 
typical hydrodynamic interactions are calculated is increased. The force coupling method 
is known to give an accurate representation of the hydrodynamic interactions when the 
distance between the beads is more than 3Aa; [2^], thus we expect to match the BD results 
for sufficiently small Ax/b. 
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FIG. 4: Normalized auto correlation functions of the Rouse-mode amplitudes, Cp{t), Eq. (j43p . 
The two longest wavelength modes {p = 1,2), an intermediate wavelength (p = 4) and a short 
wavelength (p = 8) are shown for three different grid resolutions; Ax/b = 2.58 (circles), Ax/b = 
1.29 (squares), and Ax/b = 0.86 (triangles); the solid lines are the Brownian dynamics results. 

For large grid-spacings. Ax = 2.586, the FLB diffusivity is too small, in agreement with 
previous studies [2|, while for smaller grid spacings. Ax < b, an almost exact agreement 
between FLB and BD results is found. The effect of a large grid-spacing can be understood 
from the limiting case when the grid spacing exceeds the length of the entire chain, in which 



case the polymer dynamics follows the Rouse scaling 



The recommendation 



2| that 



the grid spacing should be comparable to the mean distance between neighboring beads, 
Ax ~ (r^)^^^ = 1.606, leads to small errors in the diffusion constant, of the order of 2 — 3%. 



Rouse relaxation times 



The internal motions of the polymer coil are a more sensitive measure of the hydrodynamic 
interactions than the center-of-mass diffusion coefficient. The polymer configuration can be 
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represented by its normal (Rouse) modes Xp [32|, 



Xp = — r > r„ cos 

' ' n=0 



^ 1 

PTT 1 



P = l,2, oo, (42) 



where p denotes the mode number. The normahzed autocorrelation function, Cp{t), of the 

YJ^QdQ (p > 0), 

_ (Xp(t) . Xp(0)) 
^^(')-(Xp(0).Xp(0))' ^''^ 

decays almost exponentially 

Cp(t) = exp(-t/rp), (44) 

as can be seen in Fig.|H Tp defines the Rouse relaxation time of the p*^ mode. The correlation 
functions for the two finest grid resolution. Ax = 0.656 (not shown) and Ax = 0.866, agree 
almost perfectly with Brownian dynamics, while for the coarsest resolution Ax = 2.586 there 
are significant deviations in the long- wavelength modes. Somewhat surprisingly, the errors 
in the less resolved LB simulations diminish with increasing mode number, so that for p = 4 
the results for all grid resolutions are essentially indistinguishable. However for still shorter 
wavelengths the discrepancies increase again for the coarser grids, this time in the opposite 
direction. 

We have selected one particular time, t = 12.6to, where the p = 1 mode has decayed to 
45% of its initial value, for a more detailed comparison. For the three resolutions shown in 
Fig. m the deviation from Brownian dynamics are 11% (Ax = 2.586), 5% (Ax = 1.296), 
and < 0.5% (Ax = 0.866 and Ax = 0.656), respectively. Data in Fig. 7 of Ref. Q|, taken 
at a similar time, show deviations between FLB and BD of the order of 2%. The grid 

1 /2 

resolution in these simulations corresponds to (r^) ^ 1.1 Ax, similar to the intermediate 
resolution in our work. Ax = 1.296 or (r^)^^^ ^ 1.2Ax. Our results show slightly larger 
deviations, possibly due to the shorter chain or the softer excluded volume forces, both of 
which emphasize the short-range hydrodynamic interactions. For the system sizes we used, 
the 0{L^^) corrections to the Rouse-mode relaxation times |^ are small; simulations with 
Ax = 1.296 and L = 18.8Rg (instead of L = 9ARg) show a similar (4%) deviation from 
Brownian dynamics at t = 12.6to- 

The Rouse-mode relaxation times were calculated from the integral of the autocorrela- 
tion functions Cp{t), Eq. (jH]). The first portion of the integral was calculated by numerical 
quadrature, up to a time where the correlation function had decayed to less than 5% of its 
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Ax = 2.5Sb 



Ax = 1.2% 



Ax = 0Mb 




FIG. 5: Relaxation times of the Rouse modes, Tp, from FLB (solid circles) and BD (solid triangles). 
Results for = 10 are compared for three different grid resolutions, Ax = 2.586 (left). Ax = 1.296 
(center), and Ax = 0.866 (right); results for Ax = 0.656 are indistinguishable from Ax = 0.866. 



initial value. To minimize the effect of statistical errors on the integral, we fitted the last 
portion of the correlation function to a single exponential and calculated the long-time con- 
tribution analytically. The overall integral is insensitive to the exact value of the relaxation 
time of the fitted exponential, which was determined self-consistently from the value of Tp. 
Since the decay of Cp{t) follows a single exponential almost exactly, this procedure is quite 
precise; we used the same protocol for both FLB and BD correlation functions. 

The data in Fig. [5] show that the lattice-Boltzmann method can reproduce the whole 
Rouse spectrum when sufficiently resolved. For the two finest resolutions. Ax = 0.866 and 
Ax = 0.656 (not shown), the deviations in the long- wavelength relaxation times are less than 
1%, and are only slightly larger (~ 2%) at short wavelengths. At the intermediate resolution, 
Ax = 1.296, the long- wavelength relaxation times are well represented, with errors of 5% at 
most, but at the coarsest resolution the deviations are larger, up to 15%. Finally, in Fig. [6] 
we compare the Rouse spectrum for chains with 20 segments. The errors in the relaxation 
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Ax = 2.5Sb 



Ax=\.29b 




FIG. 6: Relaxation times of the Rouse modes, Tp, from FLB (solid circles) and BD (solid triangles). 
Results for = 20 are compared for resolutions Ax = 2.586 (left) and Ax = 1.296 (right). 

times are similar to those obtained for the shorter chains with the same grid resolution. 
IV. CONCLUSIONS 



This comparative study of lattice-Boltzmann and Brownian dynamics simulations of a 
semi-flexible polymer demonstrates that static and dynamic properties of isolated chains 
agree quantitatively, within 1 — 2%. Our paper complements recently published work 8| 
through a systematic investigation of variations in the LB model parameters. Our results 
show that hydrodynamic retardation, which is sometimes suggested to be a reason for dis- 
crepancies between LB and BD results [9^, is in fact easily controlled; the diffusion coefficient 
and Rouse spectrum are independent of Schmidt number when Scq > 10. Other parameters 
such as fluid viscosity and bead mass have little effect on the results. Somwhat disappoint- 
ingly, a higher-order interpolation of the fluid velocity field does not lead to improved agree- 
ment with Brownian dynamics. Although results with three-point interpolation converge to 
the same values as with linear interpolation, the convergence is slower, rather than faster 
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as one would have hoped. Despite the smoother interpolation of the flow fleld, the force is 
delocalized over a larger volume and this seems to reduce the accuracy, while simultaneously 
increasing the computational cost. 

The crucial parameter affecting the accuracy of a lattice-Boltzmann simulation is the res- 
olution of the polymer on the LB grid. When the mean distance between neighboring beads 
is more than twice the LB grid spacing, the agreement between FLB and BD simulations 
is essentially exact. However the computational cost of a flne grid is high, scaling as the 
resolution to the sixth^ower. Thus a reasonable practical compromise is the original sug- 
gestion (r^)^/^ ^ Ax [2]. The errors in the dynamic properties are then around 5%, which 
is sufficient for most purposes. A typical LB simulation. Ax = L296, L = dARg, a = 0.003, 
run for 8 x lO^to, requires about 70 hours of computation. Comparable Brownian dynamics 
simulations take approximately one hour. However, simulations of concentrated solutions in 
confined geometries are more favorable for lattice-Boltzmann methods j^. 
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